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Abstract. 

We investigate the relation between thermodynamic and dynamic properties of an 
associating lattice gas (ALG) model. The ALG combines a two dimensional lattice 
gas with particles interacting through a soft core potential and orientational degrees of 
freedom. From the competition between the directional attractive forces and the soft 
core potential results two liquid phases, double criticality and density anomaly. We 
study the mobility of the molecules in this model by calculating the diffusion constant 
at a constant temperature, D. We show that D has a maximum at a density p max and a 
minimum at a density p m in < Pmax- Between these densities the diffusivity differs from 
the one expected for normal liquids. We also show that in the pressure-temperature 
phase-diagram the line of extrema in diffusivity is close to the liquid-liquid critical point 
and it is inside the temperature of maximum density (TMD) line. 



§ To whom correspondence should be addressed (marcia. barbosa@ufrgs.br) 



Diffusion Anomaly in an Associating Lattice Gas Model 



2 



1. Introduction 

Water is anomalous substance in many respects. Most liquids contract upon cooling. 
This is not the case of water, a liquid where the specific volume at ambient pressure starts 
to increase when cooled below T = 4°C ilj. Besides, in a certain range of pressures, also 
exhibits an anomalous increase of compressibility and specific heat upon cooling [2j-[I]- 
Far less known are its dynamics anomalies: while for most materials diffusivity 
decreases with increasing pressure, liquid water has an opposite behavior in a large 
region of the phase diagram [5]-PB]. The increase of diffusivity of water as the 
pressure is increased is related to the competition between the local ordered tetrahedral 
structure of the first neighbours and the distortions of the structure of the first and 
second neighbours. In the region of the phase diagram where this ordered structure is 
dominant, increasing pressure implies breaking first neighbours hydrogen bonds what 
allow for interstitial second neighbours to be in a closer approach. The interactions 
are thus weakened and therefore, although the system is more dense, it has a larger 
mobility. In this sense, a good model for water and tetrahedral liquids should not only 
exhibit thermodynamic but also dynamic anomalies. In SPC/E water, the region of 
the pressure-temperature (p -T) phase diagram where the density anomaly appears is 
contained within the region of the p -T phase diagram where anomalies in the diffusivity 
are present [01 HHj . 

It was proposed a few years ago that these anomalies are related to a second critical 
point between two liquid phases, a low density liquid (LDL) and a high density liquid 
(HDL) 14J. This critical point was discovered by computer simulations. This work 
suggests that this critical point is located at the supercooled region beyond the line 
of homogeneous nucleation and thus cannot be experimentally measured. Even if this 
limitation, this hypothesis has been supported by indirect experimental results [To] ITB]. 

Water, however, is not an isolated case. There are also other examples of 
tetrahedrally bonded molecular liquids such as phosphorus [T71 [TH] and amorphous 
silica [TO] that also are good candidates for having two liquid phases. Moreover, 
other materials such as liquid metals [20] and graphite [2JJ also exhibit thermodynamic 
anomalies. Unfortunately, a closed theory giving the relation between the form of the 
interaction potential and the presence of the anomalies is still missing. 

It was observed that water has both diffusion and density anomalous behavior and 
that they are located in the same region of the p-T phase diagram. It is reasonable 
to think that the dynamics and thermodynamics are deeply related. In this case 
establishing the connection between the anomalous behavior of the diffusion constant 
and the density anomaly is a fundamental step towards understanding the source of the 
anomalies. In order to test this assumption, the dynamics of a number of models in 
which density anomaly is present were studied [221 EH] Netz et al [22] studied a system 
of molecules interacting by a purely repulsive ramp-like discretized potential consisting 
of a n number of steps of equal size. If n is small, the region on the p-T phase diagram in 
which the density anomaly is present encloses the region in which the diffusion anomaly 
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exists. If n is large, the potential becomes equivalent to a smooth ramp [23] arid the 
region p-T phase diagram in which the diffusion anomaly is present encloses the region 
in which the density anomaly exists. This behavior resembles the one expected for 
water. Unfortunately these two models do not exhibit a second critical point (or even 
a first critical point) because the interaction potential in both cases is purely repulsive. 
Therefore no connection with criticality and dynamic and thermodynamic anomalies 
could be made. 

The difficulty in including the attractive interactions in these approaches is that 
within continuous potentials attractive terms usually move the density anomaly to the 
metastable region of the p-T phase diagram. Moreover continuous potentials usually 
lead to crystallization at the region where the critical point would be expected. The 
analysis of the presence of both the two liquid phases and the critical point becomes 
indirect . 

In order to circumvent these difficulties, our model system is the a lattice gas with 
ice variables [21] which allows for a low density ordered structure [213 EE] ■ Competition 
between the filling up of the lattice and the formation of an open four-bonded 
orientational structure is naturally introduced in terms of the ice bonding variables 
and no ad hoc introduction of density or bond strength variations is needed. Our 
approach bares some resemblance to that of some continuous models |27]l28| l2§]. which, 
however, lack entropy related to hydrogen distribution on bonds. Also, the reduction of 
phase-space imposed by the lattice allows construction of the full phase diagram from 
simulations, not always possible for continuous models [22] • The associating lattice gas 
model (ALG)[2]3 |2H] is able to exhibit for a convenient set of parameters both density 
anomalies and the two liquid phases. In the framework of the ALG we address two 
questions: (i) is the presence of diffusion anomaly related to the presence of density 
anomaly? (ii) if so, what is the hierarchy between the two anomalies and the presence 
of a second critical point? We show that the two anomalies are located in the same 
region of the p-T phase diagram, close to the second critical point and that the region 
on the p-T phase diagram in which the density anomaly is present encloses the region 
in which the diffusion anomaly exists. 

2. The Model 

We consider the ALG that is a two dimensional lattice gas model on a triangular lattice 
as introduced by Henriques and Barbosa J2H] |2E]. The sites that can be empty or 
occupied. Besides the occupational variable, cr^, which assumes the value cij = if 
the site is empty, or Oi = 1 if the site is full, there are other six bond variables, t\\ 
representing the possibility of hydrogen bond formation between the neighbours sites. 
Four bonding variables are the ice bonding arms, two donor with tI° = 1 and two 
acceptor with t\ 3 = — 1. The other two arms are taken as inert, with r\ 3 = 0, and are 
always opposite to each other, as illustrated in Fig. 1. Because there are no restrictions 
on the bonding arms positions, the particles are allowed to be in eighteen different states. 
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Figure 1. Illustration of the model. Each site i can be empty or occupied and has six 
variables, rf , one for each arm. If r 4 fc = no bond is formed in spite of the configuration 
of the arm of the neighbour site, if r\ = ±1 and the arm of the neighbour site has 
Tj = =pl, a bond is formed. 



The Hamiltonian of the model has two kinds of interactions: one isotropic 
attractive, like van der Waals, and a repulsive orientational hydrogen bonding. The 
energy of the system is given by 



where -v+2u is the van der Waals interaction, -2u is the hydrogen bonds interaction, 
Oi = 0, 1 are occupation variables, rf = 0, ±1 represent the bonding arms variables as 
illustrated in Fig. 1, the summation over k is over the arms of the site we are considering 
and the summation over / is over the six neighbours arms that especific point to the 
site arm. Two neighbour sites i and j form an hydrogen bond if the product between 
their pointing arms is equal to —1, in other words, we need r*rj = — 1. In spite of each 
molecule has six neighbours, only four hydrogen bonds are allowed for site. For each 
pair of occupied sites that form a hydrogen bond, an energy -v is attributed, while for 
non-bonding pairs of occupied sites, the energy is -v+2u. 

The phase diagram of the reduced chemical potential, Jl = fi/v, versus reduced 
temperature, T = k B T/v reproduced in Fig. 2 was originally obtained by Henriques 
and Barbosa using Monte Carlo simulations for a triangular lattice with L 2 = 100 
lattice sites. At low temperatures and chemical potentials, the system is in the gas 
phase. As the the chemical is increased at constant temperature, there is a first-order 
phase transition between the gas to a low density liquid phase (LDL). As the chemical 
potential is increased even further, there is another first-order transition between the 
LDL phase a high density liquid phase (HDL). The two transition end at critical points. 

The reduced pressure, P = P/v versus reduced temperature phase diagram was 
also obtained by Henriques and Barbosa j2H| and it is illustrated in Fig. 3. In the 
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Figure 2. Phase diagram showing reduced chemical vs. reduced temperature for 
u/v = 1. The filled symbols represent the LDL-HDL coexistence line. The empty 
symbols indicate the gas-LDL coexistence line. The coexistence at zero temperature 
at ~p — —2 and ~p = 2 are exact. 
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Figure 3. Phase diagram showing reduced pressure vs. reduced temperature for 
u/v = 1. The filled symbols represent the LDL-HDL coexistence line. The empty 
symbols indicate the gas-LDL coexistence line. The coexistence at zero temperature 
at p/v = 3 and p/v = are exact. The dashed line is the TMD. 
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P — T plane, close to the LDL — HDL critical point, there is a region of temperature 
of maximum density (TMD) for a given pressure. Surprisingly the location of the 
continuous transitions, the critical points and the density anomaly are not very sensitive 
to the systems size. Some test runs were done for L 2 = 400,2500,4900. The difference 
between the critical points obtained by applying finite size scaling to those runs and the 
result for L 2 = 100 is quite small. So, for simplicity all the detailed study of the model 
properties and the full phase diagrams was undertaken for an L 2 = 100 lattice sites. 

This simple model exhibits two liquid phases, two critical points and density 
anomaly similar to the ones present in the SPC/E model for water and in other potential 
models for water. In the next section we will investigate if it also has diffusion anomaly 
and how this anomaly is related to the TMD. 



3. Diffusion 



For studying the mobility, we have performed Monte Carlo simulations of a system of n 
particles interacting as specified by the Hamiltonian Eq.Q in a triangular lattice with 
L 2 = 100 lattice sites. The procedure for computing the diffusion coefficient goes as 
follows. The system is equilibrated at a fixed chemical potential and temperature. In 
equilibrium this system has n particles. Starting from this equilibrium configuration at 
a time t = 0, each one of these n particles is allowed to move to an empty neighbour site 
randomly chosen. The move is accepted if the total energy of the system is reduced 
by the move, otherwise it is accepted with a probability exp(AE /ksT) where AE 
is the difference between the energy of the system after and before the move. After 
repeating this procedure nt times, the mean square displacement per particle at a time 
t is computed by 

(Ar(t) 2 ) = ((r(t)-r(0)) 2 ) , (2) 

where r(0) is the particle position at the initial time and r(t) is the particle position 
at a time t. In Eq. P]|, the average is taken over all particles and over different initial 
configurations. The diffusion coefficient is then obtained from the relation 

fl =l lm <^>. (3) 

*-oo At 

Since the time is measured in Monte Carlo time steps and the distance in number of 
lattice distance, a dimensionless diffusion coefficient is defined as 

D= lim ^m. ( 4) 

t^oo At V 1 

where r = r/a and a is the distance between two neighbour sites and t = tjtuc is the 
time in Monte Carlo steps. 

In order to find if our system exhibits diffusion anomalies, we have analysed how D 
varies with the number density p = n/L 2 for a fixed temperature. For each temperature 
and chemical potential, 200 samples were obtained and an average over samples was 
made. 
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Figure 5. Reduced Diffusion coefficient vs. density for T = 0.8. D increases with the 
reductions of the density as in a normal liquid. 

At low chemical potentials and T = 0.8, p is small and the mean square 
displacement has only the diffusive regime as illustrated in Fig. 4. As the chemical 
potential is increased, p also increases and the ballistic regime appears. Using the 
ballistic regime, the mean square displacement was then computed for densities ranging 
from p = 0.1 to p = 0.9 ( not illustrated in Fig. 4 for clarity). From this data together 
with Eq. (J3J), we have obtained the reduced diffusion coefficient, D, versus density, 
p, shown in Fig. 5. In this case, D increases with the decrease of the density what 
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Figure 6. Reduced Diffusion coefficient vs. density for T = 0.75. 




Figure 7. Reduced Diffusion coefficient vs. density for T = 

0.7,0.725,0.735,0.75,0.775. 

represents the behavior expected for normal liquids. Simulations for higher temperatures 
are consistent with this result. 

However, this is not the case at lower temperatures. For instance, the mean square 
displacement for T = 0.75 exhibits a peculiar behavior different from the one observed in 
Fig. 4. At the diffusive regime, for a certain density range, the slope of ((Ar) 2 ) instead 
of increasing with the decrease of density, decreases with it. Consequently the isochores 
cross each other. The reduced diffusion coefficient, D ( shown in Fig. 6) at very low 
densities decreases as p increases like in a normal liquid. However, as the density is 
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Figure 8. Reduced pressure vs. reduced temperature phase diagram showing the the 
two liquid phase, two critical points, the density anomaly and the diffusion anomaly 
regions. 
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Figure 9. Zoom of the reduced pressure vs. reduced temperature phase diagram 
showing the density anomaly and diffusion anomaly regions. 

increased D has a minimum at pDmin, and increases with the increase of density from 
PD mi „ < P < PD ma x- Increasing the density above pomaxi D decreases as in a normal 
liquid. Therefore, there is a region of densities pomax > P > PDmin where the diffusion 
coefficient is anomalous, increasing with density. This behavior is similar to the diffusion 
anomaly present in SPC/E water and ramp and shoulder models. A diffusion anomaly 
in the ALG model is observed in the range of temperatures 0.8 > T > 0.6 illustrated in 
Fig. 7. 
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The region in the p -T plane where there is an anomalous behavior in the diffusion 
is bounded by (T Dmin , Pomm) and (T Dmax , P Dmax ) and their location is shown in Fig. 8. 
The region of diffusion anomalies (To max , Pomax) and (Tn- mi-n , -Pomm) lies inside the region 
of density anomalies ( for detail see fig. 9) what differs from the behavior observed in 
SPC/E water jTOJ but coincides with the behavior shown for non smooth ramp-like 
potentials [22]- 

4. Conclusions 

In this paper we have addressed two questions: (i) is the presence of diffusion anomaly 
related to the presence of density anomaly? (ii) if so, what is the hierarchy between the 
two anomalies and the presence of a second critical point? 

For tracking the answer of these two questions we have investigated the behavior 
of the diffusion coefficient in the associating lattice gas model. This simple model is 
suitable for addressing our two questions since it exhibits two liquid phases and a line 
of density anomalies (TMD) [251 EE]- 

Using Monte Carlo simulations we computed the the mean square displacement 
with time. From the slope of this curve at the diffusive regime, the diffusion coefficient 
was derived. This procedure was implement for different temperatures and densities. 
We found that for high temperatures, particles in our model system diffuse as a normal 
liquid. At low temperatures, however, the diffusion coefficient has an interval of densities 
PD min < P < PD ma x where D increases with the increase of density. Therefore, the 
presence of a density anomaly seems to be associated with the presence of diffusion 
anomaly, confirming observations made in other models j2HHH0] and in water [TU| I12( [T^j . 
This seems to indicate that as the particles gain more energy by being close together, 
this gain facilitates the mobility. 

For addressing the second question, we have located in the p-T phase diagram 
the pressure of the density of maximum diffusivity and the pressure of the density of 
minimum diffusivity. Since for each temperature there is one PD mi „ and there is one 
PDmax > the diffusivity extrema is composed of two lines at the p-T phase-diagram. By 
comparing the TMD with the lines of diffusivity extrema, we found that the region in 
the p-T plane of diffusion anomaly is located inside the region of the density anomaly and 
that both are in the vicinity of the second critical point. This suggests that in models 
where anomalies are present, the second critical point might arises if a attractive term 
in the potential would be included. The hierarchy between the anomalies resembles 
the one observed in the purely repulsive ramp-like discretized potential [22] ■ The link 
between the two models is the presence of two competing interaction distances and 
the non smooth transition between them. The first ingredient seems to be the one 
that defines the presence of the anomalies, while the second might govern the hierarchy 
between them [22j[23j[30j. Similar behavior should be expected in other models where 
the density anomaly is also present [3*T]-|3l)] 
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